###############################################################################-
# AGEISM AND POLITICS PAPER 
# FIGURES: 1 (Generational), 2 (Personal), D6 (Both Indices)
###############################################################################-

# Clear environment and console
rm(list = ls()); cat("\014")

# Load required libraries
library(tidyverse)
library(readr)

###############################################################################-
# 1. DEFINE PATHS AND VARIABLE LABELS
###############################################################################-

# Set folder containing regression output
path <- "..../Replication_files/Output/Regressions/Tables"



# Files for each country model
files <- c("OLS_1_Grouped_cov.txt", "OLS_1_Grouped_cov_Italy.txt", "OLS_1_Grouped_cov_Korea.txt", "OLS_1_Grouped_cov_USA.txt")
countries <- c("Pooled", "Italy", "Korea", "USA")
full_paths <- file.path(path, files)

# Directory to save plots
plot_dir <- "..../Replication_files/Output/Coef_plots"

# Create save folder if it doesn't exist
if (!dir.exists(plot_dir)) dir.create(plot_dir)

# Labels for variables to use in plots
var_labels <- c(
  old_leaders_bad_index       = "Gen. Grievances Index",
  personal_grievances_index   = "Per. Grievances Index",
  anxiety_index               = "Anxiety Index",
  contact_index               = "Contact Old Index",
  ethnocentric_index          = "Ethnocentric Index",
  strong_right                = "Right-wing",
  strong_left                 = "Left-wing",
  female                      = "Female",
  age                         = "Age",
  university                  = "Bachelor or more"
)

###############################################################################-
# 2. FUNCTION TO TIDY COEFFICIENT OUTPUT
###############################################################################-

tidy_from_file <- function(filepath, country_name) {
  raw <- read_tsv(filepath, col_names = FALSE, show_col_types = FALSE)
  df <- raw[-(1:2), ] %>% filter(!if_all(everything(), is.na)) %>% mutate(row_index = row_number())
  
  coef_rows <- df %>% filter(row_index %% 2 == 1)
  se_rows   <- df %>% filter(row_index %% 2 == 0)
  
  min_len <- min(nrow(coef_rows), nrow(se_rows))
  coef_rows <- coef_rows[1:min_len, ]
  se_rows   <- se_rows[1:min_len, ]
  
  varnames <- coef_rows$X1
  labels   <- coef_rows$X2
  
  coef_data <- coef_rows %>% select(-X1, -X2, -row_index)
  se_data   <- se_rows %>% select(-X1, -X2, -row_index) %>%
    mutate(across(everything(), ~ str_remove_all(.x, "[()]")))
  
  coef_data$variable <- varnames
  coef_data$label    <- labels
  se_data$variable   <- varnames
  
  coef_long <- coef_data %>%
    pivot_longer(-c(variable, label), names_to = "col", values_to = "estimate")
  
  se_long <- se_data %>%
    pivot_longer(-variable, names_to = "col", values_to = "se")
  
  df_long <- left_join(coef_long, se_long, by = c("variable", "col")) %>%
    mutate(
      estimate   = as.numeric(str_remove_all(estimate, "\\*")),
      se         = as.numeric(se),
      col        = as.numeric(str_remove(col, "X")),
      model_col  = col,
      outcome    = case_when(
        col %in% 3:5    ~ "Ageism Explicit",
        col %in% 6:8    ~ "Ageism Thermometer",
        col %in% 9:11   ~ "Ageism Index",
        col %in% 12:14  ~ "D-Score IAT",
        TRUE            ~ NA_character_
      ),
      country = country_name
    ) %>%
    filter(!is.na(estimate)) %>%
    select(variable, label, estimate, se, outcome, model_col, country)
  
  return(df_long)
}

# Apply function to all input files
all_data <- map2_dfr(full_paths, countries, tidy_from_file)

# Add 90% and 95% confidence intervals
crit_95 <- qnorm(0.975)
crit_90 <- qnorm(0.95)

all_data <- all_data %>%
  mutate(
    ci90_low  = estimate - crit_90 * se,
    ci90_high = estimate + crit_90 * se,
    ci95_low  = estimate - crit_95 * se,
    ci95_high = estimate + crit_95 * se
  )

###############################################################################-
# 3. FUNCTION TO PLOT SELECTED MODELS
###############################################################################-

plot_and_save <- function(data, vars_to_keep, table_name) {
  
  plot_data <- all_data %>%
    filter(model_col %in% data$model_cols,
           variable %in% vars_to_keep) %>%
    mutate(
      country = factor(country, levels = c("USA", "Korea", "Italy", "Pooled")),
      variable = factor(variable,
                        levels = data$var_order,
                        labels = var_labels[data$var_order]),
      outcome = factor(outcome, levels = c("Ageism Explicit", "Ageism Thermometer", "Ageism Index", "D-Score IAT"))
    )
  
  p <- ggplot(plot_data, aes(x = estimate, y = variable, color = country, shape = country)) +
    geom_errorbarh(aes(xmin = ci95_low, xmax = ci95_high), height = 0, linewidth = 0.3, position = position_dodge(0.6)) +
    geom_errorbarh(aes(xmin = ci90_low, xmax = ci90_high), height = 0, linewidth = 1,   position = position_dodge(0.6)) +
    geom_point(position = position_dodge(0.6), size = 2.5) +
    facet_wrap(~ outcome, nrow = 1) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
    scale_color_grey(start = 0.3, end = 0.7) +
    scale_shape_manual(values = c(16, 17, 15, 18)) +
    labs(x = "Coefficient", y = NULL, color = "Country", shape = "Country") +
    guides(color = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE)) +
    theme_bw() +
    theme(
      panel.grid = element_blank(),
      strip.text = element_text(size = 10, face = "bold"),
      axis.text.x = element_text(size = 9),
      axis.text.y = element_text(size = 9),
      axis.title = element_text(size = 10),
      legend.position = "bottom",
      legend.title = element_blank()
    )
  
  ggsave(file.path(plot_dir, paste0(table_name, ".png")), plot = p,
         bg = "white", device = "png", width = 10, height = 6.5, units = "in")
}

###############################################################################-
# 4. CREATE PLOTS: FIGURES 1, 2, D6
###############################################################################-

# Figure 1: Generational Grievances
plot_and_save(
  data = list(
    model_cols = c(3, 6, 9, 12),
    var_order = c("university", "age", "female", "strong_left", "strong_right",
                  "ethnocentric_index", "contact_index", "anxiety_index", "old_leaders_bad_index")
  ),
  vars_to_keep = c("old_leaders_bad_index", "anxiety_index", "contact_index", "ethnocentric_index",
                   "strong_right", "strong_left", "female", "age", "university"),
  table_name = "Figure_1_generational"
)

# Figure 2: Personal Grievances
plot_and_save(
  data = list(
    model_cols = c(4, 7, 10, 13),
    var_order = c("university", "age", "female", "strong_left", "strong_right",
                  "ethnocentric_index", "contact_index", "anxiety_index", "personal_grievances_index")
  ),
  vars_to_keep = c("personal_grievances_index", "anxiety_index", "contact_index", "ethnocentric_index",
                   "strong_right", "strong_left", "female", "age", "university"),
  table_name = "Figure_2_personal_grievances_index"
)

# Figure D6: Both Indices
plot_and_save(
  data = list(
    model_cols = c(5, 8, 11, 14),
    var_order = c("university", "age", "female", "strong_left", "strong_right",
                  "ethnocentric_index", "contact_index", "anxiety_index",
                  "personal_grievances_index", "old_leaders_bad_index")
  ),
  vars_to_keep = c("old_leaders_bad_index", "personal_grievances_index", "anxiety_index",
                   "contact_index", "ethnocentric_index", "strong_right", "strong_left",
                   "female", "age", "university"),
  table_name = "Figure_D6_both_index"
)

